####GSEA
setwd("E:\\1.Ovarian_Treg\\1.picture\\17.GSEA")
#引用包
library(limma)
library(org.Hs.eg.db)
library(clusterProfiler)
library(enrichplot)
#读取文件,并对输入文件进行整理

rt=read.table("OV_TPM.txt", header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp), colnames(exp))
data=matrix(as.numeric(as.matrix(exp)), nrow=nrow(exp), dimnames=dimnames)
data=avereps(data)

#删掉正常样品
group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
data=data[,group==0]
data=t(data)
rownames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(data))
data=t(avereps(data))

#读取风险文件
#读取risk文件
risk=read.table("rs_tcga.txt", header=T, sep="\t", check.names=F, row.names=1)
risk$risk=ifelse(risk$riskScore>median(risk$riskScore),"high","low")
sameSample=intersect(colnames(data), row.names(risk))
data=data[,sameSample]
risk=risk[sameSample,]
#colnames(risk)[13]="risk"

Risk=risk
#高低风险比较，得到logFC
dataL=data[,row.names(Risk[Risk[,"risk"]=="low",])]
dataH=data[,row.names(Risk[Risk[,"risk"]=="high",])]
meanL=rowMeans(dataL)
meanH=rowMeans(dataH)
meanL[meanL<0.00001]=0.00001
meanH[meanH<0.00001]=0.00001
logFC=log2(meanH)-log2(meanL)
logFC=sort(logFC,decreasing=T)
genes=names(logFC)

#读入基因集文件
gmt=read.gmt("h.all.v7.5.1.symbols.gmt")

#富集分析
kk=GSEA(logFC, TERM2GENE=gmt, pvalueCutoff = 1)
kkTab=as.data.frame(kk)
kkTab=kkTab[kkTab$pvalue<0.05,]
write.table(kkTab,file="GSEA.result.txt",sep="\t",quote=F,row.names = F)
library(GseaVis)

gseaNb(object = kk,
       geneSetID = terms[5],
       newGsea = T,
       #addGene = geneList1[[1]],
       addPval = T,
       pvalX = 0.68,pvalY = 0.68,
       pCol = 'black',
       pHjust = 0,subPlot = 2,markTopgene=T,geneCol="#593327")
ggsave("5.pdf",height = 4,width = 6)

